library(tidyverse)
library(gifski)
library(ggraph)
library(here)
library(igraph)
source(here("modelFunction_rewiring.R"))

First, run the model.

# Define parameters
N = 100
edge.prob <- 0.04
burn.in = 20
burn.out = 5
pm = 0.3
ps = 0.1
pa = 0.2
add00 = c(0.5, 10)
lose01 = 0.1
add10 = 0.05
lose11 = c(0.5, 0.5)
histMultiplier = 1.2
doRemoval = TRUE
modelGraphs <- runModel(N = N, # Nodes in the network
                     edge.prob = edge.prob, 
                     burn.in = burn.in,
                     burn.out = burn.out,
                     pm = pm,
                     ps = ps, 
                     pa = pa,
                     add00 = add00, 
                     lose01 = lose01, 
                     add10 = add10,
                     lose11 = lose11, 
                     histMultiplier = histMultiplier,
                     doRemoval = doRemoval) %>%
  lapply(., function(x){
    igraph::graph_from_adjacency_matrix(x, mode = "undirected")
  })
  1. How does variation in degree emerge in this model?

Compute degree for each of the model networks.

degrees <- lapply(modelGraphs, function(x){igraph::degree(x)})
degreesDF <- do.call(cbind, degrees) %>% as.data.frame() %>% setNames(., 1:length(degrees)) %>% mutate(individual = 1:nrow(.)) %>%
  pivot_longer(cols = -individual, names_to = "slice", values_to = "degree") %>%
  mutate(slice = as.numeric(slice)) %>%
  mutate(beforeAfter = ifelse(slice <= burn.in, "before", "after"))
## Warning in (function (..., deparse.level = 1) : number of rows of result is not
## a multiple of vector length (arg 21)

Okay, so we get some degree distributions, but do we get emergent individual variation in degree that’s somewhat consistent throughout time?

# Plot everyone's degree over time
degreesDF %>%
  ggplot(aes(x = slice, y = degree, col = individual, group = individual))+
  geom_line()+
  theme_minimal()+
  geom_vline(xintercept = burn.in+1, col = "red", size = 2)

# Let's just look at a few individuals, because this plot is too hard to read
indivs <- sample(1:max(degreesDF$individual), 5)
degreesDF %>%
  filter(individual %in% indivs) %>%
  ggplot(aes(x = slice, y = degree, col = factor(individual), group = individual))+
  geom_line()+
  theme_minimal()+
  geom_vline(xintercept = burn.in+1, col = "red", size = 2)

# This is still a little hard to see... let's look at just the "before" rows.
degreesDF %>%
  filter(individual %in% indivs, beforeAfter == "before") %>%
  ggplot(aes(x = slice, y = degree, col = factor(individual), group = individual))+
  geom_line()+
  theme_minimal()

# Okay, I'm seeing some differentiation, but no clear patterns. Because of how the model is set up, the rank order tends to be maintained for a few consecutive time steps, but not over the whole span of the baseline model. This makes sense, since individuals don't have any a priori reason to have higher or lower degrees than other individuals.

Now, going to compute three network-level measures before and after the removal of an individual.

1. Network density deltas

How does edge density behave over time?

densities <- lapply(modelGraphs, function(x){
  edge_density(x)
 }) %>% 
  unlist() %>% 
  as.data.frame() %>%
  setNames(., "density") %>%
  mutate(slice = 1:length(modelGraphs))

densities %>%
  ggplot(aes(x = slice, y = density))+
  geom_line()+
  geom_smooth()+
  theme_minimal()+
  geom_vline(aes(xintercept = burn.in+1), size = 1, col = "red")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

First delta: immediately before removal to after removal/pre-rewiring

delta1_density <- igraph::edge_density(modelGraphs[[burn.in+1]])-igraph::edge_density(modelGraphs[[burn.in]])

Second delta: after removal/pre-rewiring to post-rewiring

delta2_density <- igraph::edge_density(modelGraphs[[burn.in+2]])-igraph::edge_density(modelGraphs[[burn.in+1]])

2. Network connectivity (average path length) deltas

How does connectivity behave over time?

connectivities <- lapply(modelGraphs, function(x){
  mean_distance(x)
 }) %>% 
  unlist() %>% 
  as.data.frame() %>%
  setNames(., "connectivity") %>%
  mutate(slice = 1:length(modelGraphs))

connectivities %>%
  ggplot(aes(x = slice, y = connectivity))+
  geom_line()+
  geom_smooth()+
  theme_minimal()+
  geom_vline(aes(xintercept = burn.in+1), size = 1, col = "red")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
## Warning: Removed 1 rows containing non-finite values (stat_smooth).
## Warning: Removed 1 row(s) containing missing values (geom_path).

First delta: immediately before removal to after removal/pre-rewiring

delta1_meanDistance <- igraph::mean_distance(modelGraphs[[burn.in+1]])-igraph::mean_distance(modelGraphs[[burn.in]])

Second delta: after removal/pre-rewiring to post-rewiring

delta2_meanDistance <- igraph::mean_distance(modelGraphs[[burn.in+2]])-igraph::mean_distance(modelGraphs[[burn.in+1]])

3. Network modularity deltas

How does modularity behave over time?

modularities <- lapply(modelGraphs, function(x){
  cluster_fast_greedy(x) %>%
    modularity()
 }) %>% 
  unlist() %>% 
  as.data.frame() %>%
  setNames(., "modularity") %>%
  mutate(slice = 1:length(modelGraphs))

modularities %>%
  ggplot(aes(x = slice, y = modularity))+
  geom_line()+
  geom_smooth()+
  theme_minimal()+
  geom_vline(aes(xintercept = burn.in+1), size = 1, col = "red")
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

First delta: immediately before removal to after removal/pre-rewiring

delta1_modularity <- modularity(cluster_fast_greedy(modelGraphs[[burn.in+1]]))-modularity(cluster_fast_greedy(modelGraphs[[burn.in]]))

Second delta: after removal/pre-rewiring to post-rewiring

delta2_modularity <- modularity(cluster_fast_greedy(modelGraphs[[burn.in+2]]))-modularity(cluster_fast_greedy(modelGraphs[[burn.in+1]]))

What should our measure be?

I am not at all convinced that these momentary deltas are what I should be using. Seems like we need to at least consider two time steps out, since the model relies on this. Looking at the graphs above shows me how different the longer-term trends are from the short-term fluctuations.

Let’s run the model 100 times and see if trends emerge in how to measure.

nRuns <- 100
modelRuns <- vector(mode = "list", length = nRuns)
for(i in 1:nRuns){
  modelRuns[[i]] <- runModel(N = N, # Nodes in the network
                     edge.prob = edge.prob, 
                     burn.in = burn.in,
                     burn.out = burn.out,
                     pm = pm,
                     ps = ps, 
                     pa = pa,
                     add00 = add00, 
                     lose01 = lose01, 
                     add10 = add10,
                     lose11 = lose11, 
                     histMultiplier = histMultiplier,
                     doRemoval = doRemoval) %>%
  lapply(., function(x){
    igraph::graph_from_adjacency_matrix(x, mode = "undirected")
  })
}

Compute the three network-level measures for each of the model runs.

densList <- lapply(modelRuns, function(x){
  lapply(x, function(y){
    edge_density(y)
  }) %>%
    unlist() %>%
    as.data.frame() %>%
    setNames(., "density") %>%
    mutate(slice = 1:nrow(.))
})
densDF <- data.table::rbindlist(densList, idcol = "run")

connList <- lapply(modelRuns, function(x){
  lapply(x, function(y){
    mean_distance(y)
  }) %>%
    unlist() %>%
    as.data.frame() %>%
    setNames(., "connectivity") %>%
    mutate(slice = 1:nrow(.))
})
connDF <- data.table::rbindlist(connList, idcol = "run")

moduList <- lapply(modelRuns, function(x){
  lapply(x, function(y){
    cluster_fast_greedy(y) %>%
    modularity()
  }) %>%
    unlist() %>%
    as.data.frame() %>%
    setNames(., "modularity") %>%
    mutate(slice = 1:nrow(.))
})
moduDF <- data.table::rbindlist(moduList, idcol = "run")

Make some plots and try to detect any changepoint trends

densDF %>%
  ggplot(aes(x = slice, y = density))+
  geom_line(aes(group = run, col = run), size = 0.1)+
  geom_smooth()+
  theme_minimal()+
  geom_vline(aes(xintercept = burn.in+1), size = 0.5, col = "red")+
  scale_color_viridis_c()
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

connDF %>%
  ggplot(aes(x = slice, y = connectivity))+
  geom_line(aes(group = run, col = run), size = 0.1)+
  geom_smooth()+
  theme_minimal()+
  geom_vline(aes(xintercept = burn.in+1), size = 0.5, col = "red")+
  scale_color_viridis_c()
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 100 rows containing non-finite values (stat_smooth).
## Warning: Removed 100 row(s) containing missing values (geom_path).

moduDF %>%
  ggplot(aes(x = slice, y = modularity))+
  geom_line(aes(group = run, col = run), size = 0.1)+
  geom_smooth()+
  theme_minimal()+
  geom_vline(aes(xintercept = burn.in+1), size = 0.5, col = "red")+
  scale_color_viridis_c()
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'